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Abstract 

Kinetic schemes for compressible flow of gases are constructed by exploiting the connection between Boltzmann 
equation and the Navier-Stokes equations. This connection allows us to construct a flux splitting for the Navier- 
Stokes equations based on the direction of molecular motion from which a numerical flux can be obtained. The 
naive use of such a numerical flux function in a discontinuous Galerkin (DG) discretization leads to an unstable 
scheme in the viscous dominated case. Stable schemes are constructed by adding additional terms either in 
a symmetric or non-symmetric manner which are motivated by the DG schemes for elliptic equations. The 
novelty of the present scheme is the use of kinetic fluxes to construct the stabilization terms. In the symmetric 
case, interior penalty terms have to be added for stability and the resulting schemes give optimal convergence 
rates in numerical experiments. The non-symmetric schemes lead to a cell energy /entropy inequality but exhibit 
sub-optimal convergence rates. These properties are studied by applying the schemes to a scalar convection- 
diffusion equation and the 1-D compressible Navier-Stokes equations. In the case of Navier-Stokes equations, 
entropy variables are used to construct stable schemes. 



1. Introduction 

Kinetic schemes for the compressible Navier-Stokes equations are constructed by exploiting the connection 
between the Boltzmann equation and the Navier-Stokes equations. The primary unknown in the Boltzmann 
equation is the velocity distribution function, which gives the distribution of molecules in velocity and physical 
space Q] . The macroscopic behaviour of the flow is obtained by averaging the Boltzmann equation over all the 
molecular velocities. In kinetic schemes based on finite volume or discontinuous finite elements, the intercell 
fluxes are approximated by exploiting this connection. Kinetic schemes were first developed for the solution 
of Euler equations governing inviscid compressible flows by making use of the Maxwell-Boltzmann distribution 
function, which is the equilibrium distribution for gases. While there are kinetic schemes based on a particle 
model which update the velocity distribution function, our interest in the present work is the use of kinetic 
model to construct numerical schemes at the continuum level. The kinetic model is utilized to derive a numerical 
flux function by using the upwinding principle with respect to the molecular velocities. The availability of 
such kinetic flux based schemes is useful for coupling the continuum solver with a particle-based method like 
direct simulation Monte Carlo in the case of rarefied flows where a pure continuum model would not be valid 
throughout the flow domain. The use of a DG scheme is expected to be helpful in coupling the continuum solver 
with the particle solver. Kinetic schemes at the continuum level based on finite volume discretization have been 
constructed in [21 [3] . In the approach of Mandal and Deshpande [3] , the two states on either side of a cell face 
are used to define two equilibrium distributions from which the inter-cell flux can be computed by taking the 
appropriate moments. This gives a numerical flux function which depends smoothly on the two states and can 
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be used in a method of lines approach; the resulting set of ODE are solved using any time integration scheme 
like the Runge-Kutta method or an implicit scheme. Later, the Chapman-Enskog approximation to the BGK 
model [5] of the Boltzmann equation was used by Chou and Baganoff [6] to derive a numerical flux function 
for the Navier-Stokes equations. The usual practice in finite volume methods is to use a Godunov scheme or 
approximate Riemann solver to construct the numerical flux function for the inviscid part of the flux and use 
a central difference-type scheme for the viscous part of the flux. But the use of a kinetic formulation gives a 
numerical flux function for the total flux, not just the inviscid flux. This numerical flux depends not only on 
the two states at a cell face but also on their derivatives, which are present due to the diffusion terms. 

An alternate approach is followed in the gas kinetic finite volume schemes [7J [H [5] which are based on the 
Boltzmann equation with BGK collision model. The exact solution of the BGK model is used to compute 
the time integral of the inter-cell flux in a finite volume method. Thus in this approach, the time and space 
discretizations are coupled, similar to the Lax-Wendroff scheme. The use of the BGK model automatically gives 
a scheme for the Navier-Stokes equations. Ohwada [10] has analyzed the consistency of the gas kinetic scheme, 
while a thorough analysis of the stability and accuracy of the BGK scheme for a scalar convection-diffusion 
equation has been made in [11) . 

Discontinuous Galerkin finite element methods were first proposed in [12] for the neutron transport equation. 
These methods use a discontinuous representation of the solution which allows them to compute solutions with 
steep gradients and shocks. They are useful for solving hyperbolic problems like the Euler equations of inviscid 
compressible flows and convection dominated problems like the compressible Navier-Stokes equations at high 
Reynolds numbers. The inter-cell discontinuity in the state is resolved by making use of any consistent numerical 
flux function for the convective flux. The vast developments in Godunov and approximate Riemann solvers can 
be exploited to construct these numerical flux functions. For the case of hyperbolic problems like Euler equations, 
the Runge-Kutta discontinuous Galerkin methods have been extensively developed [131 Ej and their high order 
of accuracy is demonstrated in numerical experiments. DG methods provide a rational approach to construct 
high order accurate schemes for convection-dominated problems and can be considered as a generalization of 
finite volume methods, which have been very successful in computing discontinuous solutions and convection 
dominated, high Reynolds number flows. They also have many computational advantages like a compact stencil 
which aids parallelization, and the ability to easily develop hp-adaptive algorithms. 

In the case of Navier-Stokes equations, the viscous fluxes which involve derivatives of the solution, must also 
be computed across the element boundaries. However, the derivatives are also possibly discontinuous across 
the elements due to the discontinuous numerical approximation. Averaging the derivatives to compute the 
inter-element flux leads to an unstable scheme, and is a generic problem for a DG scheme applied to any partial 
differential equation with elliptic terms. There are two approaches to construct stable schemes for the diffusion 
terms, which are sometimes refered to as the flux formulation and the primal formulation [15] . In the first 
approach, additional variables are introduced for the derivatives leading to a mixed finite element formulation. 
Some numerical flux function is required for the additional variables whose choice is dictated by the stability of 
the resulting scheme. Bassi and Rebay [TB] developed such a scheme for compressible Navier-Stokes equations 
and applied it to laminar flow problems. Cockburn and Shu [T7J developed the local discontinuous Galerkin 
(LDG) method for time dependent nonlinear convection-diffusion problems by generalizing the method of |16j . 
In the second approach, the derivatives are averaged for computing the diffusive fluxes but additional terms are 
added to achieve stability of the scheme including interior penalty terms which help to enforce continuity of 
the solution across the elements. The resulting schemes, known as interior penalty scheme, may be symmetric 
(SIPG) or non-symmetric (NIPG) with respect to the diffusive terms. Arnold [18 proposed a symmetric DG 
scheme with interior penalty for the heat equation and showed optimal error estimates in L2 norm. Baumann 
and Oden [19, 20 introduced a non-symmetric scheme for diffusion and convection-diffusion problems, which 
has been used for compressible Navier-Stokes equations in |21j . Symmetric schemes are preferable since they 
are the only ones which have adjoint consistency .15], which is crucial for obtaining optimal error estimates in 
L2 norm. A symmetric interior penalty scheme for compressible Navier-Stokes equations has been proposed 
in |22j and its optimal accuracy has been shown through numerical experiments. A symmetric space-time DG 
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scheme for Navier-Stokes equations using the appoach of [16] and diffusive Riemann solvers for the viscous terms 
has been developed in [23] and its optimal convergence rate has been demonstrated by numerical studies. The 
various DG formulations for an elliptic equation have been unified in |15j and their stability, consistency and 
adjoint consistency have been analyzed. 

The use of a kinetic scheme has the attractive feature that it provides a numerical flux function for the 
total flux which can be used in a DG scheme. The numerical flux depends on the derivatives of the solution 
which are possibly discontinuous and a naive DG discretization can lead to an unstable scheme, especially in 
the viscous dominated situation. Xu [24] has developed a DG discretization for the Navier-Stokes equations 
using the BGK schemes, which however gives sub-optimal convergence rates similar to the NIPG scheme. Liu 
and Xu [35] proposed a modified version of discontinuous Galerkin BGK scheme for Navier-Stokes equations 
which shows improved convergence rates. These BGK-based DG schemes have not yet been analyzed from the 
point of view of stability and adjoint consistent as in [15] . In the present work, we follow the interior penalty 
approach and construct stable schemes by making use of the kinetic split fluxes to add additional stabilization 
terms similar to the NIPG and SIPG schemes. The kinetic model is used to obtain a flux splitting from 
which a numerical flux function is constructed even for the diffusion terms. The development of the schemes 
is motivated by considering the scalar convection-diffusion equation for which a kinetic formulation can be 
given. Through numerical experiments on this model problem, we find the NIPG scheme gives a convergence 
rate of 0(h k ) when using Pk basis functions which is sub-optimal as compared to the usual Galerkin methods, 
while the SIPG scheme gives the optimal 0(h k+1 ) convergence rate. Similar convergence rates are observed 
through numerical experiments for the 1-D compressible Navier-Stokes equations. In the case of the Navier- 
Stokes equations, it is advantageous to use entropy variables since the equations assume a symmetric form in 
those variables H7J 12H] ■ It also facilitates the construction of entropy stable schemes which is an important 
thermodynamic constraint for compressible flows. Even the pure Galerkin scheme written in entropy variables 
satisfies the entropy condition [27] for any degree of the basis functions. For Euler equations, Barth [29, 30] has 
studied the entropy condition for discontinuous Galerkin methods and has derived a condition for a numerical 
flux function to be entropy consistent, which is a generalization of the E-flux condition [3T] to systems of 
equations. The numerical flux resulting from kinetic flux vector splitting scheme [J] satisfies this entropy 
condition. The entropy variable formulation is also advantageous since the flux jacobians are well behaved in 
the limit of incompressible flows [32] ; this is not the case with conserved variable formulation where some terms 
in the flux jacobian matrix become unbounded or assume indeterminate form. 

The rest of the paper is organized as follows. In Section ([2| we first consider the scalar convection equation 
for which a kinetic formulation is given together with a DG scheme. The energy stability of the DG scheme is 
shown for the time continuous case. In Section ^ , a kinetic formulation for the convection-diffusion equation is 
given together with DG discretization. The instability of the naive DG discretization is shown through numerical 
experiments. Stabilized versions of the schemes are constructed and their order of accuracy is studied through 
numerical experiments. In Section Q, the Navier-Stokes equations and their entropy variables are discussed. 
A novel DG scheme for the Navier-Stokes equations based on kinetic flux vector splitting is presented together 
with some numerical results which demonstrate the accuracy of the scheme. 

2. Kinetic model for convection equation 

The scalar convection equation 



is the simplest hyperbolic partial differential equation. Throughout this work, we will assume periodic bound- 
ary conditions for the above partial differential equation. The exact solution of this equation is obtained by 
convecting the initial condition at a constant speed c and the shape of the solution profile is unchanged. We 
can obtain this equation as moments of a "Maxwell-Boltzmann" distribution function 



Ut + cu x = 



(1) 




(2) 
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j — 1/2 h J + 1/2 / J+ i j + 3/2 
Figure 1: Discontinuous finite element solution 



since (<?) = u and (vg) — cu, where the angle brackets denote integration over velocity space 

<>+oo 

(■)dv 

> 

The quantity j3 > is a free parameter that can be chosen arbitrarily The "Boltzmann equation" for the 
distribution function is obtained by differentiating equation ([2| 

9t + vg x = (u t + vu x )\ — cxp \-P(v - c) 2 l 

V 7T L J 



and we use equation ([I]) to replace u t to get 

9t + vg x = (v- c)u x \ - exp [-fi(v - c) 2 ] =: Q (3) 



Note that the moment of the collision term Q is zero, (Q) = 0, so that taking moments of equation |3]) gives us 
the convection equation ([I]) . The vanishing of the moments of the collision term means that molecular collisions 
instantaneously drive the distribution function towards the local equilibrium distribution. 

2.1. Kinetic numerical flux: KFVS 

In a finite volume or DG method, we need to compute the flux across the inter-element boundaries. The 
solution is possibly discontinuous across the elements. This defines two states u + , , , u~ , , on either side of the 

1 J J+5 J+5 

cell face shown in figure (jlj, where we have used the notation 

li^ix) = lim u(x =F e) 

In the kinetic model, molecules carry information due to their streaming motion. This allows us to use the 
upwind principle to approximate the velocity distribution function at Xj + i based on the sign of the molecular 
velocity as 
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In the case of compressible gases for which kinetic theory holds, there is no free parameter in the Maxwell-Boltzmann velocity 
distribution function. A parameter like /3 is present but is related to the absolute temperature. 
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The numerical flux function is obtained by integrating the above distribution function over the velocity space 
leading to 

Fj H =F(u+ +h ,u7 H ) = (vg j+i )=F+(u+ +i ) + F-(uJ H ) 
The split fluxes are obtained by integrating over half velocity spaces, and are given by 

F ± (u) = cuA ± +uB ± , A ± = -[l±erf(s)], B ± — ± — —= exp(— s 2 ), s = cy/j3 

2 lyJ'Kp 

The flux F + is due to all molecules moving to the right (v > 0) while the flux F~ is due to all molecules moving 
to the left (v < 0); this is the upwind principle for kinetic schemes. The numerical flux function F{u~ ,u + ) is 
consistent, i.e., F(u,u) = cu and is obviously a smooth function of u. The split flux F + {u) is an increasing 
function of u while F~ (u) is a decreasing function. Hence the numerical flux F(u~ ,u + ) is a monotone flux. We 
can also write the numerical flux function as a central flux with dissipation 

F(u + ,u~) = -(cu + + cu~) + -D(u + - u~), D = cerf(s) + — = exp(-s 2 ) (4) 
2 2 y/ ir/3 

and it is easy to check that D > 0. In the limit of f3 — > oo, the split fluxes become 

i . , I cu, c > 
F + (u) = \ ' F-(u) 
10, c < 

and the numerical flux function of equation Q becomes the standard upwind flux 

F(u + ,u~) = 

We will see that the quantity D controls the rate of dissipation of energy and it is the smallest for the upwind 
scheme. In all the computations we use a value of j3 = 1. 

2.2. DG scheme 

Consider a discretization of the one dimensional domain VL by N elements Ij = [xj_i , Xj + i], j = 1, . . . ,N. 
The DG method uses the broken space of polynomials of degree k 

The DG scheme is obtained by multiplying equation ([T]) by a test function (p^ from the broken space V£ and 
integrating over element Ij 





<t> h d t u h dx - I cu h d x <t> h dx + F(u h ). j+ i(j> h (xJ k ) - F(u h ) -_i <f> h (x t ) = 

J j J ' 1 J^- 2 J 1 J 2 

The same equation would have been obtained if we had started from equation ([3| by multiplying it with a 
test function and perform the integrals over velocity space and physical space; the collision term would have 
vanished due to the integration over velocity space. Note that we have used the numerical flux function to 
resolve the discontinuity at cell boundaries, i.e., 

F(u h ) j+i =F(u+ +i ,uT +i ) 

If we take <\>h — Uh and sum over all elements, we obtain the energy equation 

Id N f N 

UhW 2 -^ cu h d x u h dx + F(u h ) i+ i luhl i+ i = (5) 

3=1 Jl * 3 = 1 
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where we have defined the inter-element jump in the solution by 



:= u h (x+ +l ) - u h (x. + i ) 



After some simple calculations, the energy equation reduces to 

JV 

2,1/-" 2 /} 



1 1 ^En! + j b ° 



i=i 



The jumps at the inter-element boundaries lead to dissipation of energy and the scheme is stable. Note that 
even if D — 0, i.e., if we use a central flux, the scheme preserves the energy for any degree of the basis functions, 
which is consistent with the exact solution of the partial differential equation. 

3. Kinetic model for convection-diffusion equation 

Consider the linear convection-diffusion equation 

Ut + cu x = [IU XX , // > (6) 
We can give a kinetic formulation for this equation using the Boltzmann equation with BGK model 

/*+«/» = — (7) 

where g is the same Maxwell-Botzmann distribution as used in the previous section. The basic idea of the BGK 
model is that the distribution function / is close to the local equilibrium distribution g, and due to molecular 
collisions, it relaxes towards g over a time interval tr which is the relaxation time. We can look for solutions to 
the BGK model by using the Chapman-Enskog expansion in the relaxation time tr 

f = fo + r R h + r\h + ■■■ (8) 

Substituting this in equation ^ and collecting terms with the same coefficient of tr, we get 

fo=9, h = -(9t + vg x ) 

as the first two terms in the expansion. The zeroth order term in the Chapman-Enskog expansion leads to the 
convection equation as seen in the previous sections. The first order term becomes 



fi = -(v- c)u x \j^ exp [-/3(v - cf 



where we have used the inviscid equation u t = —cu x to replace the time derivative term since f± should depend 
only on the lower order terms in the Chapman-Enskog expansion. The first two terms constitute the Chapman- 
Enskog distribution function 

f cc =fo + T R h = [u- tr(v - c)u x ] y|exp [-p(v - cf] (9) 
Upon taking moments of this distribution function, we obtain 

(n = u, {vn = cu-^ Ux 



G 



In order to recover the correct flux of equation ([6]), we have to choose tr = 2/3 fi so that (vf cc ) — cu — fm x . 
Since / cc is a truncated solution in the Chapman-Enskog expansion it does not satisfy equation §f§. We can 
derive an equation for / co by differentiating equation ^ as follows. 

ft* + v fx° = [«i + vu x - Tr{v - c)(u xt + VU XX )] y^ ex P [~P{ V ~ c ) 2 ] 

We now use the equation ^ to replace u t and u xt to obtain 

■ft C + v fx C = [i v - c ) u x + {n - tr{v - c) 2 )u xx - r R fi(v - c)u xxx ] y^exp [~f3(v - c) 2 ] =: Q 

Upon taking moments of this equation, the right hand side vanishes, i.e., (Q) = 0, and we obtain the convection- 
diffusion equation 

3.1. Flux vector splitting 

The split fluxes for the viscous problem are obtained by integrating the Chapman-Enskog distribution 
function f cc over the half velocity spaces, — oo < v < and < v < +oo. The expressions for the kinetic split 
fluxes can be written as 

F ± = (cu - ^iu x )A ± + uB ± 

The numerical flux function F(u + , u~) = F + (u + )-\-F~ (u~) can be written as a sum of convective and dissipative 
fluxes 

F = F c {u+,u-) + F d (ut,u-) 
where the convective flux F c is identical to equation Q, while the dissipative flux is given by 

F d (u+,u-) = -[1U+A+ - fxu~A- 

Note that the diffusive flux can also be written in split form as 

F d (u+, u-) = F+(u+) + F-(u-), F±{u x ) = -nUvA* 

In the case of the convection-diffusion equation the viscous fluxes depend only on the derivatives of the solution, 
but for the Navier-Stokes equations, the viscous fluxes will depend on the solution also. 

3.2. Test case 

We will consider the convection-diffusion equation ([6| in the domain (— 1,+1) together with the initial 
condition 

u(x, 0) = — sin(7rx) (10) 
and periodic boundary conditions. The exact solution to this problem is given by 

u(x, t) = — exp(—fiir 2 t) sin(7r(a; — ct)) (11) 

In all the tests, we set the convection speed to be c = 1. We define two test cases based on different values 
of the viscosity coefficient [i as given in table ([lj with tf being the final time. Test case 1 corresponds to a 
situation in which convection dominates diffusion, while in Test case 2, diffusion dominates convection effects. 
The time integration is performed using the 3-stage, strong stability preserving Runge-Kutta scheme of Shu 
and Osher [53]. The time-step is chosen based on a CFL condition of the form At = CFL • m&x(h/\c\, h 2 /n). 
In the case of convection equation solved with a DG scheme, the CFL number depends on the degree of the 
polynomial basis functions [34] . 
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Test Case 


1 


2 


M 


0.001 


1 


*/ 


30 


0.5 



Table 1: Parameters for the test case based on convection-diffusion equation K\ 



3.3. DG scheme 

The DG discretization of equation (|6| is obtained by multiplying it by a test function (f>h and integrating 
over any element Ij 

4> h d t u h dx- / cu h d x (j) h dx + /i / d x u h d x <fi h dx + F(u h ) 1+ i<f> h (x+ 1 )-F(u h ) j _i(j) h (x~ 1 )=0 (12) 
.// .// 2 

In the above equation, the inter-element flux is approximated by the numerical flux which is defined as 

F(u h ) = F c (u+,u^) + F d (d x ul,d x ul) 

and it includes both convective and diffusive fluxes. Substituting (f)^ = Uh and summing over all elements we 
obtain 

Id N f N f N 

2 eh 7 ""' 1 " 2 ~~ ^ / CUhdxUhdx + / (9 x u h ) 2 dx + y^ y F(u h ) j+ i = 

Following similar steps as in the case of linear convection equation, we obtain the following energy equation 

* N N r N 

^^H 5 2 + i / (d x u h ) 2 dx + ^2F d (d x u h ) j+i K] i+ i =0 (13) 



1 d II II 2 

2dt IM 2' 



The second and third terms are dissipative since they lead to a decrease of the energy. But we do not know 
the sign of the last term, which can be either positive or negative. We apply this DG scheme to the two test 
problems using Pi basis functions and the results are shown in figure ([2f. When /i = 0.001 the results look 
correct but this is not so when fi — 1 which corresponds to a viscous dominated case. The scheme is not able 
to dissipate energy at the correct rate. We see that a naive discretization of convection-diffusion equation by 
the discontinuous Galerkin method leads to unstable schemes when viscosity is dominant. 

3.4- DG scheme with stabilization 

Motivated by the energy analysis, we propose the following stabilized DG scheme 



") h d t u h dx - / cu h d x (f) h dx + n / d x Uhd x ^hdx + F(uh) i+ i(l>h(xJ,i)-F(uh) i _i<f>h(x. 1 ) 

Jij Jij J 5 3 5 (14) 

- F+(d x <t>+) J+h {u h \ jH - F^id^l)^ M^x = 

Note that the last two terms in the above equation are constructed from the kinetic split viscous fluxes. Choosing 
4>h = Uh and summing up over all elements gives us 



1 A N N 

1 d II 112 

2dt KI1 _ , 

3=1 J=l"^ 



N N 

^^H J+ i+pE / (d x u h fdx = (15) 



which shows that the scheme is stable. The choice of the additional stabilization terms in equation (14) was 
just enough to eliminate the problematic term in the energy equation (13). The stabilized scheme is applied to 
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Figure 2: Convection-diffusion problem using equation i ] 1 2[ i and Pi basis functions, (a) Test case 1, (b) Test case 2 

convection-diffusion problem with ^ = 1 and the results are shown in figure [3j The stabilized scheme is able to 
dissipate the energy better than the unstabilized one, but it exhibits a significant amount of phase error for Pi 
basis functions; the energy dissipation is still seen to be incorrect with Pi basis functions. However, the use of 
Pi basis functions is seen to give excellent results without any phase error. In the next sections, we see that 
adding a penalty term gives correct results even with P\ basis functions. 



3.5. Cell energy analysis 

In the previous sections, we performed a global energy analysis; here we make a cell energy analysis. If we 
multiply the convection-diffusion equation by u and integrate over element Ij, we obtain 



dt 1 



1 



-cu 



J- 2 



+ [-ftuuj^i = -/i / (u x ) z dx < 

J 2 



(16) 



We will now try to mimic this property for the DG scheme for which we will have to define consistent numerical 
fluxes for the second and third terms on the left hand side of the above equation. Substituting (ph = Uh m 
equation (14 1, we obtain 

d 



df 



Wfn 



u h( x j + i) - u hOj_i)J + F ( u h) j+ iu h (x+ + ,) - F{uh)j_iUh(x,_i) 
- F+{d x ul) j+ i - Fj(d x u^) j _i Inkjet = -H I {d x u h ) 2 dx 



The numerical flux F(uh)j + i consists of the convective and viscous fluxes. Define the average value 

{{w/Jj+i := \{uh{x+ + i) +u(xj + i)) 
and convective and diffusive numerical fluxes 

F c (u h ) j+ i = F c {u h ) J+ i lu h } j+ i - - {{u 2 h }} j+ i 
F d (u h ) J+ i = Ff(d x u£) j+ iu h (xJ +i )+F4(d x Uk) j+ iu h (x+ + ±) 
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Figure 3: Results for Test case 2 using equation \14\ : (a) Pi and (b) P2 basis functions 



The last two equations provide consistent numerical fluxes for ^cu and —[iuu x , respectively. Then we obtain 
the cell energy equation for the DG scheme as 



which is consistent with the exact energy equation (16). The inter-element jumps in the solution add additional 
dissipation in the numerical solution. 

3.6. DG scheme with interior penalty 

The above examples have shown the importance of energy stability in constructing good DG scheme for the 
convection-diffusion equation. The schemes were written for each element but it is more instructive to write 
a weak formulation for the whole computational domain since it helps to understand some structure of the 



equations. By summing up the equation ( 14) for each element, we write the DG scheme for convection-diffusion 
equation as a weak statement: Find Uh{t) £ Vt such that for all 4>h £ V£, the following equation is satisfied 



4> h d t u h dx- I cu h d x 4> h dx + fj, / d x u h d x <p h dx + V] F c (u h ) j+ i l<p h j ., i 
Jn Jn ^ 2 

+ ^2F d {d x u h ), J+ i_ l<ph\ J+ i +eJ2 F d(d x <l>h) j+ i + Yj S '^ Uh h+k 



(17) 



1 = 



where the last term is an additional term that has been added, called the interior penalty term, and is of the 
form [TBI [H] 

S h (u) = C ip ^luj 

where C- lp is a non-dimensional positive number. Two classes of schemes are obtained based on whether e — — 1 
or +1. Define the bilinear form consisting of all diffusive terms 



a e (vh,w h ) = y^ J F d (d x v h ) i+ } 



i+i 



eJ2 F d(d x w h ) J+ i 



i+i 
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1. If we choose e = — 1 in ( 17) and omit the interior penalty term, then we obtain the scheme of equation ( 14 1. 
The scheme is non-symmetric since a—i(vh,Wh) ^ o,-i(wh,Vh)- As we have seen, this scheme satisfies a 
cell energy inequality. The interior penalty terms are necessary to make the scheme stable since 



N 

a-i(v h ,v h )=H / (d x v h ) 2 dx + Q 



JV 

lp h 



and the right handside is a norm on the space of discontinuous functions only in the presence of the penalty 
term for any C lp > 0. Without the penalty term (C- lp = 0), the right hand side is only a semi-norm and 
the scheme is only weakly stable 15J. In this case, the discretization of diffusive terms is similar to that 
of Baumann and Odcn 20 for elliptic equation, who also show sub-optimal convergence in the case k > 2 
only. There is no proof of convergence for k = 1 case which combined with only weak stability explains 
the poor performance of the scheme shown in figure (|3^,). With the interior penalty terms, the global 
energy equation is 

1 a 1 N N r N 



2dt' lUh 



.7=1 ~ .7=1 Jl i 3=1 



which shows that interior penalty adds additional stabilization to the scheme which is important when 
viscosity is dominant. The scheme with interior penalty is known as NIPG and is stable for any k > 1. 
2. If e = +1, the discretization of the diffusive terms becomes symmetric, i.e., ax{vh,Wh) — Q>i(wh,Vh), an d 
the scheme is known as SIPG scheme. However, the SIPG scheme does not satisy a cell energy inequality, 
though in the case of elliptic equation it can be shown to be globally stable |15j . This scheme is stable 
only with the addition of a penalty term 5^, with a sufficiently large value of C lp . The determination 
of a good value of C lv is an open problem; in practice a value of C; p between 5 to 10 seems to yield 
satisfactory results. A larger value of C\ p reduces the allowable time step for explicit time integration 
schemes, requiring the necessity of implicit time integration schemes. 

The jump terms in the above DG scheme do not destroy the consistency of the scheme; if we substitute the 



exact solution of the convection-diffusion equation which is continuous in equation (17), it exactly satisfies the 
equation since the jump terms vanish. 



3. 7. Numerical results using interior penalty scheme 

In this section, we present results for the two test cases, using the interior penalty scheme given by equa- 
tion ( fl7| ). We first present results using 20 cells. The schemes NIPG and SIPG are tested with a penalty 
parameter C\ v = 10. Figure Q shows the results for Test case 1 which is convection dominated situation. 
Figure ([5| shows results for Test case 2 which is a diffusion dominated problem. We see that the addition of 
the interior penalty leads to better solutions with P\ basis functions, as compared to figure |3]). 

In order to compute the convergence rates of the various schemes, wc perform a grid refinement study. 
Errors in the solution and its derivative are computed in L 2 norm and shown in tables ([2j9| for P\ and P^ basis 
functions. We observe the following trends: 

• SIPG gives more accurate results than NIPG. This can be seen by comparing the error values for the two 
schemes given in the tables. 

• NIPG converges at sub-optimal rates; both the solution and its derivative converge at 0(h k ). 

• SIPG converges at optimal rates; the solution converges at 0{h k+1 ) and the derivative converges at 0(h k ). 
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Figure 5: Results for Test case 2 using interior penalty scheme 
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Table 2: Test case 1: NIPG, Pi basis 
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Table 3: Test case 1: NIPG, P 2 basis 



Actually both the schemes seem to exhibit slightly higher convergence rates on Test case 1 which is convection 
dominated. In Test case 2 which is viscosity dominated, the convergence rates are very close to the above 
mentioned values. These convergence rates are consistent with what is observed in the case of DG schemes 
applied to elliptic equation [15]. The NIPG scheme for elliptic problems is numerically found to show optimal 
convergence rate for odd degree polynomials but we have not observed this phenomenon for the convection- 
diffusion equation with the current NIPG scheme. 

3.8. Test case 3 

Torrilhon and Xu |11| have analyzed the BGK finite volume scheme for scalar convection-diffusion equation 
and find that the BGK scheme exhibits only first order accuracy in the asymptotic limit, though on coarse grids 
it shows third order accuracy. They find a non-monotone convergence of the error with grid size. In order to 
study the behaviour of the error convergence, we apply the present DG scheme on a test case similar to the one 
given in [11, with /3 = 1, fi = 0.005, c = 1, using Pj basis functions. The grid is refined starting with 50 cells 
and increased by 50 cells in each refinement step. The initial condition is taken to be 

8 16 
u(x, 0) = 4 + - sin(7re/2) H sin(37ra;/2) 

7T 37T 

for which an exact solution is available |llj for the case of periodic boundary conditions using Fourier approach. 
The convergence of the error in the solution as a function of the number of cells is shown in figure (16a) 
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Table 4: Test case 1: SIPG, Pi basis 
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Table 5: Test case 1: SIPG, P 2 basis 
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Tabic 6: Test case 2: NIPG, Pi basis 
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Tabic 7: Test case 2: NIPG, P 2 basis 
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Table 8: Test case 2: SIPG, Pi basis 
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Table 9: Test case 2: SIPG, P 2 basis 
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for Pi basis functions; this plot must be compared with figure (7) in [TT]. It is clear that both the NIPG and 
SIPG schemes converge monotonically with grid refinement. Asymptotically, the NIPG scheme converges at 
first order while the SIPG scheme converges at second order. Surprisingly, the NIPG scheme exhibits second 
order convergence on coarse grids, a superconvergence phenomenon that is also observed in the BGK finite 
volume scheme The convergence of the error for P2 basis functions is shown in figure ([6J3) which shows 

second and third order convergence for NIPG and SIPG schemes, respectively. On coarser grids, the NIPG 
scheme again shows optimal convergence rate. These convergence rates are consistent with what is observed in 
the previous two test cases, and also in the 1-D Navier-Stokes equations given in the later sections. 



4. Navier-Stokes equations 

The Navier-Stokes equations which represent conservation laws for mass, momentum and energy for a fluid 
can be written as a system of partial differential equations 

d t U + 8iFi(U) + diGi{U, VU) = 

where U is the vector of conserved variables corresponding to density of mass, momentum, and energy, while 
Fi, i = 1,2,3 are the inviscid or convective fluxes and Gi are the viscous fluxes. The conserved variables are 
given by 

U = [ p pui pu 2 pu 3 pe ] 

where p is the density, u — (ui, 1(2,113) is the fluid velocity and e is the energy per unit mass, which can be 
written in terms of the internal energy e and the kinetic energy as e = e + \u 2 . An equation of state e = e(p,p) 
is required to close the system of equations^ where p is the pressure; for a polytropic ideal gas, the equation 
of state is given by e = > where 7 = j¥- is the ratio of specific heats at constant pressure and constant 



1G 



volume |35j . The flux vectors are given by 



F, 



pu % 







pSn + pu. k ui 




-Til 


p5 i2 + pUiUi 


j Gi — 


-Til 


p6 a + puiU 3 




-m 


(pe+p)ui _ 







(18) 



The shear stress tensor r is related to the strain rate by Newton's constitutive law while the heat flux is related 
to the gradient of absolute temperature T by Fourier's law 



djUi) 



-fJ,(d k u k )Sij , 



-nd,T 



(19) 



The coefficient of dynamic viscosity p and the coefficient of heat conduction n are related through the Prandtl 
number Pr = We will use the notation F = [F%, F 2 , F 3 ] and G = [G\, G2G3]. Moreover, for any vector 

ii e I 3 , we denote F ■ n = FiUi and G ■ n = Giiii for the fluxes in the direction n. Since the numerical 
experiments in this paper are conducted for the 1-D Navier-Stokes equations, the corresponding equations are 
given in |Appcndix A| 



4-.1. Euler equations, entropy variables and symmetric hyperbolic form 
The system of equations containing only the convective terms 

d t U + diFi = 



(20) 



are called the Euler equations, and are hyperbolic in nature. This means that for any vector n = (m, n<2, 713), 
the matrices 

A(U, n) = Ai(U)ni, MU) = F!(U) 

have all real eigenvalues and a complete set of eigenvectors. The Euler equations admit an entropy function 
r)(U) with entropy fluxes 0i(U), i = 1,2,3 such thai[^] 

e' i {u) = v '{u)F' i {u) 

This implies that for smooth solutions, an additional conservation law is satisfied, i.e., d t r\ + d$i = 0, while for 
discontinuous solutions the inequality d t r\ + diOi < is satisfied in the weak sense. The entropy and entropy 
fluxes can be related to the physical entropy s = \n{p/ p 1 ) + sq by 



rj = 



PS 

7-1' 



pSUi 



1 



(21) 



Actually any function of the form 77 = ph{s) where h is convex is an entropy for the Euler equations, but the 
above linear choice is the only one which is consistent with the second law of thermodynamics in the presence of 
heat conduction |27j . Due to the mass conservation equation, the constant term Sq can be ignored and moreover 
it can be shown that rj(U) is a strictly convex function. We can then define the entropy variables by 



V T = rj'(U) 



(22) 



which can be inverted in a unique manner to obtain U = U(V). Let us define the dual variables by a Legendre 
transform 

r,*(V) = V ■ U(V) - T)(U(V)), 9*(V) = V ■ Fi(u(y)) - e % {u{v)) 



2 For a scalar function of a vector like rj(U) we consider v'(U) to be a row vector. 
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Differentiating these expressions it is easy to check that 77* (V) = U T and 9* (V) = Fi(U (V)) T . It then follows 
that the matrices 17' (V) = rf (V) and Fl(U(V))U'(V) = 6*" (V) are symmetric. Moreover the matrix 

U'(V)=r,"(U(V))- 1 



is positive-definite. Making a change of variables from U to V in equation (20), we get 

A d t V + AAV = 

where we have defined 

Ao = U'(V) = rj*"(V), Ai = Ado = 0*"(V) 
the matrix A is symmetric positive definite, while the matrices Aj are symmetric. 

4-2. Kinetic description of Euler equations 

The Boltzmann equation pQ gives a molecular statistical description of gas dynamics in the form of the 
velocity distribution function f(v,x,t) which can be written as 

dtf + v t d t f = Q(f,f) 

where the right hand side represents the effect of intermolecular collisions. The meaning of the distribution 
function is that the quantity f(v, x, t)dxdv gives the average number of particles found in the spatial volume 
dx around the point x at time t and having velocities in the range dv around the velocity v. Thus / gives the 
number density of molecules in phase space and by a suitable scaling, it can be considered as providing the 
mass density. The macrosopic fluid motion is the average effect of all molecular motions, so that the conserved 
variables are obtained by integrating the distribution function over the velocity space 



U= / 1 p(v)f(v,x,t)dv = (iPf) 
7r3 

where 

ip(v) = [ 1 vi v 2 v 3 i|w| 2 ] 

are the microscopic collisional invariants. Since collisions do not destroy mass, momentum and energy, we have 
(ipQ) = 0. In fact, it can be shown pQ that any function <f>(v) of the molecular velocities for which (<fiQ) = 
must be of the form <f>(v) = a i4>i- 

A gas in thermodynamic equilibrium is characterized by the fact that there are no spatial or temporal 
variations in the distribution function, i.e., Q(f, f) = 0. It can be shown that the unique solution g = f of this 
equation is of the form 

g = cxp(V ■ ip) 

for some constants KgK 5 . This is the usual Maxwell-Boltzmann distribution function and can be written in 
the more familiar form 



3/2 



2RT 

The Maxwell distribution can also be characterized as the solution of the following minimization problem (36) 

min {H[f] : (V/) = U} , H[f] = (f In /) 
and the quantities V appear as the Lagrange multipliers to enforce the constraint. If we define 

V*(V) - (/) = (exp(V-VO) 

then 

rf'(V) = (i/> T exp(y- VO) = U T (23) 
and hence the Lagrange multipliers V are precisely the entropy variables corresponding to the macroscopic state 



U. The modifications in the case of polyatomic gases are briefly discussed in Appendix B 
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Figure 7: Definition of edge normal and trace values 



4-3. From Maxwell distribution to Euler equations 

The Maxwell distribution was obtained under the condition of thermodynamic equilibrium corresponding to 
a gas which has uniform properties in space and time. The appropriate moments of the distribution function 
yield the conserved variables and the inviscid fluxes, i.e., U = {^{v)g) and F t = (vi^(v)g). Due to this 
correspondence, it is common to use the collision-less Boltzmann equation 

d t g + Vidig = 

to construct numerical schemes for the Euler equations. When the gas is not uniform, then the assumption is 
that locally, it is still in thermodynamic equilibrium, so that one has a local Maxwell distribution 

g(v,x,t) = exp(V(x,t) ■ -0) 

defined at every point of space. This distribution however does not satisfy the collision-less Boltzmann equation. 
We can derive an equation for g by differentiating the above equation 



-A 1 A i + Vil 



d t g + Vidig = g [d t V + v l d l V] ■ ip = g 
Taking moments on both sides, we obtain 

(iP(d t g + vAg)) = -(ip <» ^A^AAV + (v^ ® ipg)diV (24) 



Now from equation ( 23 ) 

n*"{V) = {i>®i>g) = U'{V) = A Q 

Moreover, since Fi = (vitpg), we get F[(V) — (Vitp^ipg) — Ai, which shows that the right hand side of equation 
(24) vanishes, and we obtain the Euler equations. Due to this reason the right hand side is ignored and we 
use the collision-less Boltzmann equation with Maxwell-Boltzmann distribution function to derive numerical 
schemes for Euler equations. 

4-4- DG scheme for Euler equations 

Consider a triangulation Th = {K} of the domain by polygonal, non-overlapping elements. We assume 
that the edges of the triangulation are oriented, so that each edge has a unit normal vector n with trace values 
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of the discontinuous functions on the edge defined as in figure (j7|). Let us also define the jump and average 
operators as 

H = + -(•)"= l} = \[(-) + + (•)"] 

respectively. Note that we do not involve the normal vector n in the definition of these operators since it is used 
in the definition of the numerical fluxes. The DG scheme for an interior element K is given by multiplying the 
Euler equation by a test function Wh and doing some integration by parts to obtain 

/ W h ■ d t U(V h )dx - f Fi(V h ) ■ diWhdx + f H C (V+, V h ~,n) ■ W+da = (25) 

JK JK JdK 

The integral on the element boundary dK is to be interpreted as the sum of the edge integrals, i.e., J K (-)da = 
Seedif J (-)d(T. Note that the flux across the element boundaries is approximated by a numerical flux function. 
The numerical flux is given by the concept of upwinding at the kinetic level; since molecules carry mass, 
momentum and energy, the flux across any edge is determined as a contribution due to molecules crossing the 
edge in both directions, which leads to the following numerical flux function 

H C (V+, V^,n) = ((v ■ n) + ^g{V+)) + ((v ■ n)"^(V^)) (26) 

where we have used the notation ^ 

(v-n)* =-[(v-n)± \v-n\] 

Note that each of the integrals in the two terms of the numerical flux are effectively taken over half the velocity 
space, i.e, v ■ n > and v ■ n < 0, respectively. The kinetic flux can also be written as 



H c (V h + , V h -,n) = \ ((v ■ n)iPg(V+)) + ±((v n)^g{V^)) + \(\v r# [g(V+) - g(V^)]) 



and using 



ds 

we can write the numerical flux as 



(K) - 9(V h ~) = I 1 ^-g(V( S ))d S = [ g(V(s))^ ■ lV h j ds, V(s) = sV+ + (1 - s)V^ 



^c{V^V h -,n) = \ [F{V+)+F{V-)]-n+ l -J^ (\v ■ n\ij) ® ipg{V(s))) ■ {V h } ds 

This is precisely the kinetic mean- value (KMV) flux of Barth [371 ED]- Hence for the Euler equations, the 
KMV flux is identical to the KFVS flux, which can be written in terms of explicit formulae but involving error 



functions and exponentials [3]. The expressions for the one dimensional split fluxes are given in the Appendix 


4-4-1- Entropy stability of DG scheme 



To study entropy stability of the scheme, we follow Barth and substitute Wh = Vh in equation (25) to get 
f d tV (V h )dx = - [ [Q(V+, V h ~ , n) + D(V+,V h -,n)] da 

JK JdK 



where 

e(v+,v-,n) = ivj-n c (v + 1 v-,n) - 
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is a consistent and conservative numerical flux for the entropy flux 8 ■ n, and 



D(V+,V-,n) = \ (lVj-n c (V+,V-,n) - [0* • n|) = \ \ \V\ ■ [U C {V+, V~,n) - F(V(s)) ■ n) ds 



For entropy stability, we have to show that D is positive; a sufficient condition is given by [141 130] 

[VI- [H c (V+,V-,n)-F(V(s))-n] > 0, Vs G [0, 1] (27) 

which is a generalization to systems of equations, of Oshcr's E-flux condition [IJT] for scalar conservations laws. 
Barth [30] has shown that the KMV flux and hence also the KFVS flux, satisfies the E-flux condition; thus the 
KFVS-based DG scheme for the Euler equations satisfies the entropy condition for any degree of basis functions. 

4-. 5. Chapman- Enskog distribution function 

Consider the Boltzmann equation with the BGK model for the collision term 

d t f + v i d i f=^^ (28) 

If we write the solution as a series in the relaxation time as in equation ^ , we obtain the following first two 
terms of the solution 

fo = 9, h = ~{dt9 + v t d t g) 

The distribution function g depends on the macroscopic variables p, it, T; we substitute the inviscid equations 
governing these quantities into the expression for f\ to obtain the following Chapman-Enskog velocity distribu- 
tion function [S] 

where the relaxation time is given by tr = ii/p and = Vi — Ui is the peculiar velocity of the molecules. It can 
be checked that U = (^/ ce ) and Fi + Gi — (vi4>f cc ), i.e., the moments of the Chapman-Enskog distribution 
yield the conserved variables and fluxes for the Navier-Stokes equations. Since f cc is a truncated solution in 



r = g 



(29) 



the Chapman-Enskog expansion, it does not satisfy the Boltzmann equation (28); we can derive an equation for 
/ ce as dtf ce + Vidif cc — Q and it can be shown through explicit but lengthy computations that (tpQ) = 0. The 
Chapman-Enskog distribution is derived for a monatic gas for which Prandtl number is one; however written 



in the form of equation ( 29 ) , the effect of Prandtl number is accounted in the computation of the heat flux q 



and we take the above distribution function to be valid for polyatomic gases also [51 [35] . 

4-6. DG scheme for Navier-Stokes equations 

In terms of the entropy variables, the Navier-Stokes equations can be written as 

A d t V + AAV - diiKijdjV) = 0, d = -K^djV 

and the matrices Kij are symmetric positive semi-definite [39| 127] . Taking dot product with entropy variables, 
we obtain 

dtv + OA + di(V ■ d) = -(diV) J KijdjV < (30) 

which is essentially the entropy condition. In order to mimic the above property, it is useful to write the DG 
scheme for Navier-Stokes equations using entropy variables. For an interior element K, the DG scheme which 



is the analoque of the scheme given by equation ( 14 1 is 



W h ■ d t U(V h )dx - / [F t (V h ) + Gi(V h , W0] • diW h dx 

K JK ^ 31 ^ 

/ H(V+,VV+, V h -, VV h - ra) • W+da - [ H+(V+, VW+,n) ■ [V h j da = 

JdK JdK 
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where n is taken to be the outward normal to the element K. The numerical flux function % which consists of 
both inviscid and viscous fluxes is given by kinetic splitting as 

H(V + , W+ , V- , W _ , n) = H C {V+, V~,n) + H d (V + , W+ V~ , VV~ , n) 



Here W c is the KFVS numerical flux for Euler equations as given by equation ( 26 1 while Hd is the KFVS 
numerical flux for the viscous fluxes obtained from Chapman-Enskog distribution. This flux can be written as 

H d (V + ,VV + ,V~ ,VV~ ,n) = Hj(V + , W+,n) (V~,VV~,n) 



Note that equation (31 ) is a non-symmetric DG scheme for the NS equation at the element level. 
4-. 6.1. Entropy stability 



In order to study entropy stability of the scheme given by equation (31) we substitute Wh — Vh to obtain 
d tV (V h )dx + f [Q(V+, V h ~,n) + D(V+, V-,n)] da 



K 



OK 



+ / [HZ(Vjt, VV+,n) ■ V,; + Ha(V^,VV h -,n) ■ V+] da - f G t (V h , VV h ) ■ diVdx - 

JdK JK 

The first two integrals are common with the Euler equations, while the last two terms are the contributions 
from NS equations. We observe that 

• V h -+Hj (V h -,VV h -,n) ■ V+ 

is a consistent and conservative numerical flux for the viscous entropy flux V ■ (GiUi), while from the Navier- 
Stokes equations, we have 

Gi{y h ,VV h ) ■ diV h = -(diVhfKijdjVh < 



which leads to a cell entropy inequality for the DG scheme which mimics equation (30). The non-symmetric 
nature of the scheme (31) is required to obtain the above cell entropy inequality. This type of local or cell 
entropy condition does not seem to be valid in the case of the symmetric version of the scheme. 

4-6.2. DG scheme for NS equation 

We now state the DG scheme for NS equation in a global formulation including penalty terms which is 
analogous to the scheme given in equation (17). Let T% denote the set of all interior edges/faces of the finite 
element grid and let T denote the boundary faces. The time continuous DG scheme for Navier-Stokes equations 
is given by 



W h ■ d t U{V h )dx 



[Fi{V h ) + Gi{V h , Wh)] • diW h dx 



nV h + , V h -,VV h -,n) ■ {Wh} da 



+e f H d (V+, VW+, V-, VW h - n) ■ \V h \ da + f S h {V h ) ■ {Wh} da + N r {V h , W h ) = 
Jrv Jr-r 



(32) 



Choosing e = — 1 yields the NIPG scheme while e = +1 yields the SIPG scheme. The interior penalty term is 
of a similar form as in the scalar problem and given by 



Sh{V h ) = C ip —{U{V h )\ 
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where we use the coefficient of kinematic viscosity in order to have dimensional consistency. The boundary 
conditions are imposed through the boundary terms Np in a weak form |22j and given by 

N T (V h ,W h ) =jH(TV+,TVV+,TV^,TVV^,n) ■ W+da 

+ J n d (TV+,TVW+,TV,;,TVWf;,n) ■ [V+ - TV+]da + J S r (V h ) ■ W+da 

The operator T modifies the boundary trace of the solution to apply the boundary conditions; for example, in 
the case of a no-slip boundary 

= rv h - = [v+ ,0,0,0, v+] T , r w+ = r w fc - = w+ 

For an isothermal wall, is modified by using the specified wall temperature. For an adiabatic wall, VV^ 
is modified so that the heat flux is zero while at a farfield boundary with conditions Voo, we take TV } + = , 
TV^ - = Voo and rVV^~ = rVV^ = VV^~. The boundary penalty term Sr is similar to the interior penalty and 
is of the form 

S r ( Vh ) = C ip ^[U(V+)-U(TV+)} 

4-7. Test case: Order of accuracy 

The NS equations do not have 1-D analytical solutions which could be used to compute the convergence 
rate of the numerical schemes. Hence we construct exact solutions by the method of manufactured solutions. 
We assume the following form for the exact solution 

p{x) = 1 + -cos(27ra), u{x) = 10x 2 (l - xf sin(2vrx), T(x) = 1 + 2x 2 (l - xf (33) 

which solves the following stationary 1-D NS equations with source terms 

(pu)x = fi> (p + pu 2 )x - t x = f 2 , {{pe+p)u) x - (tu) x + q x = f 3 (34) 

together with the boundary conditions u(0) = u(l) = and q(0) = q(l) = 0, which are the adiabatic boundary 
conditions corresponding to zero heat flux. The source terms /i, /2, fa are computed by substituting the exact 
solutions into the above equations. The viscosity coefficient is p = 0.01 or p = 1 which corresponds to convection 
dominated and viscous dominated cases, while the ratio of specific heats 7 = 7/5 corresponding to polyatomic 
gas like air. Computations are performed on a sequence of four uniform grids with N = 40, 80, 160 and 320 
elements using Pi , P2 and P3 basis functions. The L2 norm of the error in velocity and temperature and their 
derivatives are computed from the exact solution by performing a high order quadrature. The convergence 
of the error in velocity and temperature and their derivatives are plotted in figure ^ and ^ for the case of 
p = 0.01 and p = 1 respectively. The short lines are drawn with slopes indicated beside them. From these 
figures we can observe that for the NIPG scheme, the solution and its derivative converge approximately at the 
rate of 0(h k ). For the SIPG scheme, the solution converges at the rate 0(h k+1 ) while the derivative converges 
at the rate 0(h k ). The SIPG scheme is slightly better in terms of the error in the solution while both schemes 
yield nearly similar solution for the derivatives. 

4-8. Test case: 1-D shock structure 

We compute the shock structure using the Navier-Stokes equations. While it is known that the NS equations 
are not valid inside the shock, we use this as a verification test case. The exact solution of the NS equations 
for a stationary shock are obtained using the method of [JHj which is implemented as a matlab code in [3]. 
This solution is not truly exact since it uses an ODE solver, and hence we are not able to compute errors in 
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Figure 10: Density for shock structure problem using entropy variables 



the DG solution. The parameters defining the problem are: Mach number ahead of the shock is Mi = 1.5, 
7 = 5/3, Pr = 2/3 while the viscosity law is given by fi = /xi(T/Ti) ' 8 where the subscript "1" denotes pre-shock 
conditions and /Zi = 0.0005. The time discretization is made with backward Euler implicit scheme and the time 
step is chosen based on the convective speeds as At — CFL • max (j^| +e ) with CFL = 5. Figure |Io| shows the 
results for the density using NIPG/SIPG schemes and P1/P2 basis functions on grids with N = 100, h = 1/400 
and TV = 200, h — 1/800 elements, while figures (11) and (12 1 show the corresponding results for shear stress 
and heat flux. It is clear that the numerical solutions are able to compute the shock structure quite accurately. 
It is more important to resolve the variation of shear stress and heat flux since these exhibit an extremum inside 
the shock. In the case of Pi basis functions, the solution derivative is constant inside each cell and so we plot 
the constant value of shear stress and heat flux at the center of each cell with a symbol. It is clear from these 
plots that both the schemes give accurate solutions and are able to resolve the shock structure well, even for the 
shear stress and heat flux. The DG solution compares very well with the analytical solution of the NS equation. 
With P2 basis functions, we see that the DG solution is able to capture the variation of shear stress r and heat 
flux q inside each cell very accurately; the accuracy is seen to improve with grid refinement. 



5. Summary and conclusions 

We have proposed a DG scheme for Navier-Stokes equations which makes use of kinetic flux vector splitting 
to construct a numerical flux function for both the convective and diffusive fluxes. These schemes are motivated 
by applying them to scalar convection-diffusion problem and studying their accuracy. Energy/entropy stable 
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Figure 11: Shear stress for shock structure problem using entropy variables 
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Figure 12: Heat flux for shock structure problem using using entropy variables 
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schemes are obtained by adding additional stabilization terms which are constructed from the kinetic numerical 
flux and interior penalty terms. The resulting schemes which can be either symmetric or non-symmetric in 
the discretization of diffusive terms, exhibit good accuracy properties. The symmetric schemes show optimal 
convergence rates in numerical tests on the scalar convection-diffusion equation. These good properties are 
observed also for Navier-Stokes equations based on entropy variables. Future work with these schemes will be 
focused on applying them for two dimensional viscous flow problems to study if the observed convergence rates 
carry over to higher dimensions, and trying to couple them with particle methods for rarefied flow situations. 



Appendix A. 1-D Navier-Stokes equations 

In conservation form the 1-D Navier-Stokes equations can be written as 

8U OF dG 
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and the shear stress and heat flux are given by 
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The coefficient of heat conduction is given hyn 
Prandtl number. 



Pr 
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(7-1) p 



- where R is the gas constant and Pr is the 



Appendix B. Polyatomic gases 

The Maxwell-Boltzmann distribution function discussed above corresponds to a monatic gas for which the 
ratio of specific heats is 7 = | . Real gases like air which are mostly composed of diatomic molecules N2 and O2 
have a value of 7 = 1.4 and effectively behave like a diatomic gas. The internal energy of such gases is different 
since they have contributions from rotational degrees of freedom. In order to get the correct value of 7 and hence 
the internal energy, additional internal energy variables can be introduced. In the approach of Deshpande |28j , 
the collisional invariant corresponding to energy is modified to \\v\ 2 + I while in the BGK schemes [9], it is 
modified as \\v\ 2 + £ 2 , with £ 2 = £ 2 + . . . + the value of K depending on 7. A third approach is to use 
\\v\ 2 + I s as introduced by Perthame |36j : this has the advantage of preserving the form of the Boltzmann 
entropy function H = f In / [351 [3U] and will be adopted in the present work. Note that the kinetic split fluxes 
resulting from any of the above approaches are identical since there is no splitting with respect to the additional 
degrees of freedom. With the above choice of the collisional invariants, the Maxwell-Boltzmann distribution 
function is given by |30j 



j(p,u,T;v,I) 



a(i?T) 1 /(7-i) 
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1 dl, 



6 = 
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7-1 



The moments of the above distribution function yield the correct expressions for the energy and energy flux for a 
polyatomic gas while the mass and momentum equations are unaffected. The entropy variables for a polytropic 



gas corresponding to the convex entropy (21 ) are given by 



V = 



s \u\ u 1 

^T~2£T' RI" ~RT 
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Any constant terms in the entropy variables, e.g., in the definition of the physical entropy s, can be ignored 
since they do not change the entropy stability properties of the NS equations or of the DG schemes since only 
the jump in the entropy variables or their derivatives appear in the equations. 



Appendix C. KFVS fluxes for NS equation in 1-D 



The inviscid kinetic split fluxes are given by 



puA ± + pB ± 
(j> + pu 2 )A ± + puB ± 
(pe + p)uA ± + (pe + p/2)B 



where 




A ± 




1 



exp(-s 2 ) 



The viscous kinetic split fluxes are given by 



- [t + Ipuq] f3B ± 
-tA ± + pqB ± 
{~uT + q)A ± - [Ir+lPuq] B 



where r and q are the shear stress and heat flux. 
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